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1  Statement  of  the  problem  studied 

The  detection  and  discrimination  of  unexploded  ordnances  (UXO)  is  a  challenging  task  which  re¬ 
quires  sophisticated  forward  models  and  efficient  inversion  algorithms.  Due  to  the  large  number 
of  unknowns  (target’s  size,  composition,  position,  orientation,  and  number),  the  latter  often  do  not 
solve  the  entire  non-linear  problem  at  once.  Instead,  specific  but  optimized  techniques  are  typically 
developed  for  the  inversion  of  a  subset  of  parameters,  and  are  used  in  conjunction  with  one  another. 

Although  the  determination  of  the  target’s  size  and  composition  usually  requires  advances  math¬ 
ematical  models  -  which  have  been  the  subject  of  some  of  our  previous  work  —  the  determination  of 
the  position  and  orientation  is  more  tolerant  to  some  approximations  and  hence,  can  be  performed 
in  a  simpler  manner.  Typically,  an  optimization  approach  is  often  used,  whereby  the  measured  mag¬ 
netic  field  is  matched  to  the  one  produced  by  a  forward  model,  in  which  the  unknown  coefficients  are 
solved  for  by  a  matrix  inversion.  This  approach,  however,  can  produce  an  ill-posed  inverse  scattering 
problem  when  applied  to  the  inversion  of  the  location  and  orientation  of  the  object. 

Recently,  an  alternative  method  to  accomplish  this  inversion  without  having  to  solve  an  ill-posed 
problem  has  been  proposed.  The  method  assumes  that  the  secondary  field  induced  by  the  buried  tar¬ 
get  is  akin  to  that  of  three  orthogonal  dipoles,  whose  locations  and  moments  are  obtained  analytically 
from  the  estimation  of  the  magnetic  vector  potential  and  the  electric  scalar  potential,  in  addition  to 
the  measured  magnetic  field.  These  potentials  are  computed  via  the  intermediate  step  of  introducing 
surface  magnetic  dipoles  over  a  fictitious  boundary,  whose  amplitudes  are  determined  by  matching 
the  magnetic  fields,  and  relating  these  dipoles  to  the  potential  via  well  known  relations.  Once  the 
three  quantities  are  known,  the  distance  vector,  the  eigenvalues  of  the  polarizability  tensor,  and  the 
Euler  angles  are  obtained  by  a  simple  and  well-behaved  matrix  inversion. 

Following  this  approach,  we  investigate  the  robustness  of  this  new  method  to  various  amounts 
of  noise  in  the  measured  magnetic  field  as  well  as  the  influence  of  various  numerical  parameters  on 
the  accuracy  of  the  estimated  position  and  dipole  moment. 

The  major  drawback  of  the  previous  approach,  however,  is  that  it  cannot  be  generalized  to  mul¬ 
tiple  targets  (at  the  time  of  this  writing).  Yet,  multiple  targets  are  often  encountered  in  the  field, 
with  UXO  buried  alongside  other  UXO  or  alongside  clutter  items.  It  is  therefore  important  to  de¬ 
vise  an  analysis  method  flexible  enough  to  address  these  more  challenging  situations  as  well.  With 
this  purpose  in  mind,  we  propose  and  formulate  an  inversion  algorithm  based  on  an  iterative  Gauss- 
Newton  method,  for  which  the  number  of  targets  is  in  principle  unlimited.  Initial  implementation 
and  validation  of  the  method  are  presented. 
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2  Summary  of  the  most  important  results 

2.1  (H,A,y/)  method  to  estimate  an  object  location,  orientation,  and  dipole  mo¬ 
ment 

2.1.1  Background  and  procedure 

Our  past  efforts  have  focused  on  developing  some  of  the  most  up-to-date  forward  and  inverse  algo¬ 
rithms  for  the  modeling  of  UXO  in  the  electromagnetic  induction  (EMI)  regime,  in  both  the  frequency 
domain  and  the  time  domain.  These  include  the  use  of  the  spheroidal  system  for  the  general  mod¬ 
eling  of  axially  symmetric  UXO  [1—3]  (a  more  general  model  than  the  dipole  model  [4, 5]).  the  use 
of  analytical  metrics  for  the  characterization  of  UXO  [6],  the  use  of  learning  algorithms  for  their 
discrimination  [7],  and  the  recent  generalization  of  our  modeling  capabilities  to  the  ellipsoidal  coor¬ 
dinate  system  [8]. 

Although  the  detection  of  UXO  is  relatively  simple,  their  discrimination  is  a  vastly  more  com¬ 
plex  task.  Treating  as  unknowns  the  location,  orientation,  and  composition  of  the  object  produces 
highly  non-linear  problems  which  are  both  challenging  and  time-consuming  to  solve.  Effectively,  the 
forward  algorithms  mentioned  previously  become  drastically  more  efficient  if  some  unknowns  can 
be  determined  by  add-on  methods.  Based  on  a  dipole  approximation,  an  analytical  method  has  been 
recently  proposed  for  the  inversion  of  the  location  and  strength  of  the  dipoles,  effectively  providing 
information  on  the  location  and  orientation  of  the  target  in  the  ground  [9],  The  method  is  based  on 
the  near- field  approximation  of  a  magnetic  dipole,  for  which  the  field  and  the  two  potentials  are  given 
by  [10]: 


H  = 


AkR 


4ttR3 


3 R(R  ■  m) 

R |2 


m  (1  -i kR)-kzR  x(R  x-m) 


A  =  /j.0m  Jl(l-ikR)elkR 


y/  = 


4; tR3 
R  ■  m 


4nR3 


(1  -i  kR)e 


ikR 


(la) 

(lb) 

(lc) 


where  R  =  r-r',f  and  r'  being  the  position  of  the  observation  and  of  the  source,  respectively.  Note 
that  we  have  not  made  any  assumptions  as  to  which  term  could  potentially  be  neglected,  but  certainly 
some  terms  can  be,  given  that  we  work  in  the  limit  of  k  — *  0.  Typically,  the  k2  term  in  Eq.  (la)  can  be 
neglected,  as  well  as  (1  -iAR)exp(iAff),  which  leaves  a  problem  independent  of  k,  i.e.  independent  of 
frequency. 

These  equations  can  be  inverted  to  yield  R  and  m  as  function  of  H,  A,  and  i//,  and  yield: 


-  _  2 Hy/  +  (H  xA/p q) 

\H\ 2 


(2a) 

(2b) 


where 

ikR 

G(R)=  —  (1-ikR)*— .  (3) 

4  n  4  n 

The  whole  purpose  is  therefore  to  estimate  A  and  tfr  knowing  only  the  magnetic  field  H  over  a 
measurement  grid. 
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2.1.2  Determining  R 


Following  [9],  we  use  a  fictitious  surface  of  magnetic  dipoles  defined  by 

N  N  N 

Mr')  =  *E  PXnK(r)  +  y  E  PnFyn{f')  +  z  E  P^(r'). 

71  =  1  11  =  1  11=1 


(4) 


The  functions  Fx,y,z  are  basis  functions  whereas  the  unknown  coefficients  Pn’y,z  are  determined  by 
matching  the  measured  (known)  vectorial  magnetic  field  over  a  reference  grid.  Starting  from  Eq.  (21) 
in  [9],  we  write  the  magnetic  field  as 


N  N  N 

H(f)  =  E  Pnfx(r  -  0+  E  Pynfy(r  -  r'n)+  E  P*nfzir  -  r'n) 

71  =  1  71=1  11  =  1 


(5) 


where 


fn(r 


r'n)=  ^  dS  (3( 


r  —  rre>[(r  —  rre)’7?]  —  lr  —  2" 


Ir-r'J5 


(6) 


and  where  f]  =  x,y  ,z  and  Sn  refers  to  the  discretization  of  the  fictitious  surface.  From  Eq.  (5),  the  x 
component  of  the  magnetic  field  is  given  by: 


N  N 

x-H(r)  =  Hx(r)=  E  fX(F  ~Ki)+  Y  PnX- fy(f  ~r'n) 

71  =  1  71  =  1 

N 

+  YPnX-fZ(r-f'n).  (7) 

71  =  1 

Evaluating  this  equation  over  P  observation  points  (typically  over  the  measurement  grid): 


Hx(f  i) 

x-fx{r1-r'l) 

x-rin-r'j 

••• 

px 

ri 

Hx(r2) 

x-fx{f2-r\) 

x-fx(f2-r'2) 

•••  x-fx(r2~r'N) 

px 

r2 

Hx(fP ) 

x-fx{fP-f\) 

x-fx(rP-r'2) 

•••  x-fx(rp~r'N) 

px 

r  ivJ 

x-fy{ri-r'1) 

x-fy(r\-r’2 ) 

...  x-f_y(fi-f'N) 

>?' 

■ 

x-  fy(r2-r'1) 

x  •  PXr  2  -  r'2) 

...  x-fy(r2-r’N) 

H 

x  ■  fy{rp  -  r\) 

x-  fy(rp  -  r'2) 

...  x-fy(rP-r'N) 

PN 

x- fz{ri-r'2) 

...  x-  fz(ri  —  r'N) 

P{ 

+ 

x  •  fz(r2  -  f'l) 

x-fz(f2-r'2) 

...  x-  fz(r2  —  r'N) 

Pz 

r2 

x-fz(rP  -r'j) 

x-  fz(rP  -r'2) 

...  x-  fz(rP  ~r'N) 

Pz 

:N\ 

This  equation  is  rewritten  in  matrix  form  as 


Hx=FxxPx+FyxPy+Fzx-Pz. 

(9a) 

Similar  equations  can  be  written  for  Hy  and  Hz : 

=  Fxy  Px  +Fyy  Py  +Fzy  Pz, 

(9b) 

Hz=Fxz-Px+Fyz-Py+Fzz-Pz. 

(9c) 

-3- 


Delpsi,  LLC 


Prepared  for:  Army  Research  Office  (contract  #W911NF-08-l-0388) 


Real  dipole  moments 

m  =(1,0,2) 

Real  dipole  position 

(x,y,z)=(  10  cm, -2  cm, -45  cm) 

Fictitious  surface 

z'  =  -30  cm 

Measurement  grid 

8x8  points  over  [-50,50]  cmx[-50,50]  cm 

Fictitious  grid 

8x8  points  over  [-50,50]  cmx[-50,50]  cm 

Table  1:  Definition  of  the  reference  parameters. 


The  system  can  eventually  be  gathered  in  a  matrix  form  as 


Hx 

P_XX 

Pyx 

P_zx 

px 

Hy 

= 

Fjxy 

F_yy 

P_zy 

py 

(10) 

Hz\ 

P xz 

Pyz 

Pzz. 

Pz 

or 

H  =  PP, 

from  which  P  can  be  obtained,  typically  as  a  least  square  solution: 

p  =  [PT -Py1  -{pT -h). 


(id 


(12) 


Once  these  coefficients  are  obtained,  the  dipole  moments  over  the  fictitious  surface  are  obtained 
from  Eq.  (4),  from  which  (H,A,y/)  are  obtained  via  Eqs.  (1),  from  which  the  position  of  the  dipole  is 
obtained  via  Eq.  (2a). 


For  the  numerical  illustration  of  the  method,  the  parameters  of  Table  1  are  used,  unless  spec¬ 
ified  otherwise.  The  variations  in  the  retrieved  positions  due  to  other  parameters  is  illustrated 
subsequently.  Note  that  in  the  subsequent  results,  we  plot  r'  =  r—R.  In  addition,  the  captions  of  the 
figures  list  the  values  of  the  root  mean  square  errors  (rmse)  of  Ft, 


PH „  = 


PH  =  (PHx,pHy,PHz), 


Tj  =  x,y,z 


(13a) 

(13b) 


and  similarly  for  A  and  for  the  one  value  of  V.  The  root  mean  square  errors  give  some  sense  as 
to  what  components  are  well  or  poorly  estimated,  and  their  influence  on  the  position  estimation. 
Typically,  the  magnetic  field  rmse  should  be  very  low  (less  than  10-4  since  it  is  the  field  used  to 
derive  the  unknown  coefficients*),  whereas  the  vector  and  scalar  magnetic  potentials  should  have  an 
rmse  of  10~2  or  less. 

The  first  result,  which  can  be  used  as  a  reference,  is  shown  in  Fig.  1,  and  illustrates  the  fact 
that  the  retrieval  of  R  is  not  trivial.  In  particular,  important  oscillations  are  seen  in  the  retrieved 
components.  Nonetheless,  the  average  values  are  relatively  reasonable,  and  the  variances  are  be¬ 
tween  1  cm  and  2  cm,  which  would  be  acceptable  in  a  field  measurement.  The  influence  of  other 
parameters  such  as  the  position  of  the  fictitious  surface,  the  size  of  the  measurement  grid,  the  size 
of  the  fictitious  grid,  or  the  use  of  some  components  of  the  magnetic  field  only  are  all  illustrated  in 
Figs.  2-6. 


*Note  that  this  order  of  magnitude  does  not  necessarily  hold  when  only  few  components  of  the  magnetic  field  are  used. 
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Figure  1:  Retrieved  position  for  the  parameters  of  Table  1.  The  amplitudes  of  the  oscillations  are  quite  large 
but  the  variances  are  reasonable  (compared  to  the  field  precision  required).  The  average  values  are  relatively 
close  to  the  expected  ones,  p#  =(1.3, 1.1,1. 7)xl0-9,  pA  =(1.4, 2.5,2. 7)xl0-2,  py  =5xl0-3. 


x’avg=6.6  [em],  ox=6.392  [em]  y’aVg=-1.1  [cm],  ay =4.861  [cm]  z’gvg=-45.5  [cm],  oz=4.971  [cm] 


(a)  z'  =  -20  cm:  Overall  the  estimation  of  the  location  is  worse,  with  in¬ 
creased  oscillations  in  all  the  components  (reflected  also  by  larger  variances). 
pH  =(0.94, 0.60, 1.06)xl0“12,  pA  =(1.5,2.9,3.0)x  10-2,  py  =3.3xl0~2. 


x’  =-0.6  [cm],  a  =27.932  [cm]  y’  =-0.1  [cm],  a  =24.680  [cm]  z’  =-29.6  [cm],  a  =5.740  [cm] 

avg  L  J  x  L  J  avg  1  J  y  1  J  avg  z  1  J 


(b)  z'  =  -10  cm:  The  estimates  become  significantly  worse.  It  is  probably  due 
to  the  fact  that  being  so  far  from  the  real  source,  the  fictitious  dipoles  do  not 
capture  the  original  field  very  well.  Increasing  the  size  of  the  fictitious  surface 
and  the  number  of  fictitious  dipoles  seems  to  improve  the  results  (see  Fig.  5). 
pH  =(5.19,3.76,5.04)xl0“15,  pA  =(0.59, 1.24, 0.66)x  10“\  py  =  1.79x10“  ^ 


=10.1  [cm],  ct  =2.699  [cm] 


=-44.6  [cm],  o  =2.995  [cm] 


fl/WW^ 


0  20  40 

Grid  point 


(c)  z'  =  -44  cm:  The  estimate  is  relatively  good,  except  in  y' .  The  good  re¬ 
sults  are  probably  due  to  the  proximity  of  the  fictitious  source  to  the  real 
source,  thus  capturing  most  of  the  radiated  field.  pH  =(3.73,2.27,4.36)x  10-5, 
pA  =(3.6,2.6,4.7)xl0-2,  py  =2.16e-04. 


Figure  2:  Varying  the  position  of  the  fictitious  surface  z' . 
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Tm 

W*1 

HI 

0  100  200 

Grid  point 

Grid  point 


Grid  point 


(a)  16  x  16  points  over  [-50,50]  cmx[-50,50]  cm:  The  results  are  relatively  sim¬ 
ilar  to  the  reference  Fig.  1.  pn  =(4.39,3.83,5.05)xl0-4,  pa  =(1.7,2.7,3.3)x  10-2, 
pv  =4x  10-3. 


=7.7  [cm],  a  =3.572  [cm] 


=-2.4  [cm],  a  =4.665  [cm] 


=-44.8  [cm],  a  =2.694  [cm] 


10 
I  0 
$  -10 


-20, 


0  100  200 
Grid  point 


0  100  200 
Grid  point 


(b)  16  x  16  points  over  [-100,100]  cmx  [-100,100]  cm:  Increasing  the  size 
of  the  measurement  grid  does  not  necessarily  improves  the  results.  Most 
likely  the  fictitious  grid  has  to  be  adjusted  as  well.  pH  =(0.9,0.8,1.4)x  10-3, 
pA  =(0.9,1.2,1.6)xl0-2,  pv  =lx  10-3. 


*’avg=7.7  [°ml’  o,=3.503  [cm] 


y’avg=-2.4  [cm],  cry =4.578  [cm] 


z’avg=-44.9  [cm],  az=2.655  [cm] 


-30i - 

-60o - s 

0  10C 

Grid  point 


(c)  32x32  points  over  [-100,100]  cmx[-100,100]  cm:  Same  conclusion  as  the  previ¬ 
ous  case  (note  that  a  difference  inversion  scheme  has  been  used,  because  of  com¬ 
putation  time).  pH  =(1-0, 0.8,1. 3)x  10-3,  pa  =(0.9,1.3,1.6)xl0-2,  py  =lxl0-3. 


Figure  3:  Varying  the  size  of  the  measurement  grid. 
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x’a  =1 0-1  [cm],  0x=1 .1 79  [cm]  y’avg=-1 .6  [cm],  oy= 1 .437  [cm]  z’avg=-45.8  [cm],  az= 1 .692  [cm] 


Grid  point  Grid  point  Grid  point 


(a)  16  x  16  points  over  [-50,50]  cmx  [-50,50]  cm:  Somewhat  similar  results 
to  the  reference  Fig.  1.  pn  =(5.19,5.48,7.61)xl0-10,  pA  =(1.4,2.7,2.9)x  10-2, 

pv  =4x  10-3. 


x’avg=13.9  [cm],  ox=6.572  [cm]  y’aVg=-1.5  [cm],  oy=9.129  [cm]  z’avg=-42.3  [cm],  oz=5.218  [cm] 


(b)  16  x  16  points  over  [-100,100]  cmx[-100,100]  cm:  The  results  are  degraded 
(essentially  seen  in  the  larger  variances).  Probably  due  to  a  well  over-determined 
system.  pH  =(4.04, 4.66, 6.16)x  10“10,  pA  =(3.0,7.8,8.4)x  10"2,  pv  =9xl0"3. 


Figure  4:  Varying  the  size  of  the  fictitious  grid. 


X'avg=7‘5  [cm]’  ax=23-373  tcm] 


y’avg=-0.5  [cm],  oy=23.144  [cm]  z’avg=-28.3  [cm],  oz=3.613  [cm] 


Figure  5:  Increased  fictitious  grid  when  further  away  from  the  real  source:  the  results  are  improved  compared 
to  Fig.  2(b).  pH  =(8.33,5.57,8.98) xlO-16,  pA  =(0.50, 1.65, 0.95)xl0“1,  pv  ^^xlO”1. 
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x’avg=3.6  [cm],  ox=1 1.877  [cm]  y’avg=-0.9  [cm],  cy=9.583  [cm]  z’avg=-40.4  [cm],  oz= 3.222  [cm] 


(a)  Using  only  (Hx,Hy)  (thus  a  sub-system  of  Eq.  (10)).  The  results  are  signifi¬ 
cantly  degraded,  pn  =(0,0,8-15)xl0-2,  p a  =(1.8,2.8,2.9)xl0-2,  py  =6.9xl0-2. 

x’avg=11-5  [cm],  ox=7.791  [cm]  y’avg=-1.2  [cm],  oy=9.055  [cm]  z’avg=-39.9  [cm],  oz=3.752  [cm] 


(b)  Using  only  Hz  (thus  a  sub-system  of  Eq.  (10)).  The  results  are  better  than  us¬ 
ing  only  (Hx,Hy),  suggesting  that  standard  GEM-3  measurements  might  be  rea¬ 
sonably  usable.  p}j  =(4.5,2.1,0)xl0-2,  pa  =(2.6,8.9,6.9)xl0-2,  py  =3.1xl0_1. 

Figure  6:  Using  only  some  components  of  the  magnetic  field. 
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2.1.3  Determining  m 

The  average  dipole  moments  are  determined  from  Eq.  (2b)  using  the  estimated  positions  and  poten¬ 
tials.  The  results  are  summarized  in  Table  2,  along  with  those  for  the  average  R,  for  the  same  case 
as  in  the  previous  section.  The  table  illustrates  the  influence  of  various  parameters  on  the  accuracy 
of  the  inverted  results,  in  particular  the  size  of  the  measurement  grid,  the  side  of  the  fictitious  grid, 
and  the  number  of  components  of  the  magnetic  field  used  in  the  inversion. 

2.1.4  Using  GA  to  determine  R  and  m 

The  encouraging  results  of  the  previous  section  prove  that  Eqs.  (2)  and  the  use  of  a  fictitious  bound¬ 
ary  for  unknown  3D  dipoles  provide  a  good  way  to  estimate  the  location  and  moments  of  a  dipole 
knowing  only  the  magnetic  field  over  an  observation  grid.  In  real  life,  however,  results  cannot  be 
compared  against  a  known  truth  like  previously,  where  the  exact  position  of  the  dipole  is  known. 
Hence,  another  metric  needs  to  be  chosen. 

In  what  follows,  we  use  the  standard  deviation  in  the  estimate  of  R  ( i.e .  axyz  in  the  previous 
results)  as  a  metric,  supposing  that  the  smaller  it  is,  the  more  accurate  the  estimate.  This  is  some¬ 
what  intuitive  because  we  expect  the  mean  to  be  more  trustworthy  if  the  data  are  more  concentrated 
around  it. 1  This  metric,  which  exact  formulation  is  examined  subsequently,  provides  a  measure  of 
goodness  so  that  we  can  minimize  it  using  some  stochastic  methods.  We  have  implemented  here  a 
Genetic  Algorithm  (GA)  for  that  purpose,  using  5  unknown  parameters: 

(i)  The  starting  coordinate  xmin  of  the  observation  grid  in  x,  within  the  bounds  -60  cm  and  -40  cm. 
We  suppose  that  the  observation  and  fictitious  grid  overlap  in  the  ( xy )  plane  and  are  only 
different  by  their  z  position. 

(ii)  The  starting  coordinate  ymin  of  the  observation  grid  in  y,  with  the  same  assumptions  and  same 
bounds.  In  addition,  we  suppose  the  grid  to  be  1  mx  1  m  so  that  the  first  two  parameters  totally 
define  the  observation  and  fictitious  grids  in  the  (xy)  plane  as  [xmin,xmin  +  1  m],  [ymin,ymin  + 
lm], 

(iii)  The  vertical  position  of  the  fictitious  grid,  within  -50  cm  and  -10  cm. 

(iv)  The  number  of  measurement  points  N„  on  the  observation  grid  (same  in  x  and  y). 

(a)  The  number  of  dipoles  N([  on  the  fictitious  grid  (same  in  x  and  y). 

The  other  GA  parameters  are  a  population  size  of  100,  a  probability  of  crossover  of  60%,  a  probability 
of  mutation  of  5%,  and  20  iterations.  The  choice  of  the  metric  is  likely  to  change  the  estimates  of  R 
and  m,  so  that  we  examine  this  effect  in  Table  3.  For  example,  if  the  metric  is  chosen  to  be  only  oz, 
the  fitness  of  the  reference  result  is  1.572  xl0~2,  as  shown  in  Fig.  1. 

The  results  again  show  the  influence  of  various  parameters  on  the  accuracy  of  the  results.  In 
particular,  the  measured  dataset  (all  components  of  the  magnetic  field  or  only  its  z  component)  are 
seen  to  have  an  important  impact  on  result  accuracy.  It  is  expected  that  the  new  generation  of 
EMI  sensors  will  be  able  to  measured  full  vectorial  magnetic  fields  with  multiple  primary  incidences, 
which  should  improve  field-based  inversion  results  significantly. 


'  Note  that  this  assumes  that  the  estimator  is  unbiased,  which  we  expect  to  be  since  Eq.  (2a)  is  an  exact  relation. 
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Table  2:  Retrieved  average  positions  and  moments  for  the  cases  in  Section  2.1.2. 
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Table  3:  GA  results  for  various  metrics. 

'  Since  the  fictitious  surface  is  below  the  estimated  source  depth,  this  case  is  re-evaluated  by  forcing  the  fic¬ 
titious  surface  at  -46.2  cm,  to  yield  the  new  parameters:  R  =  (11.93,-1.98,-45.7)  cm,  m  -  (0.84,0.002,2.02), 
fitness=1.8xl0-2. 

*  We  proceed  in  the  same  way  as  in  to  obtain:  R  =  (11.6,-2.17,-45.5)  cm,  m  =  (0.83,0.002,1.92), 
fitness=2.54x  10-2. 
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2.2  Detection  of  two  targets  using  a  Newton  based  method 
2.2.1  Motivation 

Although  the  previous  (H,A,y/)  method  has  been  shown  to  yield  reasonable  position  and  dipole  mo¬ 
ment  estimates,  we  have  found  that  its  major  drawback  is  that  it  is  not  straightforward  to  generalize 
to  a  multi-target  configuration,  whereby  multiple  dipoles  (position  and  dipole  moments)  need  to  be 
inverted  for.  Yet,  multiple  targets  are  often  encountered  in  the  field,  with  UXO  buried  alongside 
other  UXO  or  alongside  clutter  items. 

In  the  remainder  of  this  report,  we  propose  and  formulate  an  alternative  method  to  address  this 
issue.  The  method  is  based  on  an  iterative  Gauss-Newton  implementation  which,  in  principle,  does 
not  limit  the  number  of  unknowns  and  hence,  can  potentially  invert  for  multiple  positions  and  dipole 
moments  simultaneously.  As  we  show  subsequently,  the  Jacobian  matrices  required  by  the  method 
take  analytical  forms  within  a  dipole  approximation,  which  we  therefore  follow  in  order  to  enhance 
the  efficacy  of  the  method.  Initial  results  on  two-target  inversion  are  shown  to  be  encouraging. 

We  assume  a  bistatic  configuration  were  two  targets  are  buried  at  positions  r'x  and  r’2  and  re¬ 
spond  with  a  single  dipole  moment  each,  denoted  by  mi  and  m2,  respectively.  The  purpose  is  to  use 
the  measured  components  of  the  magnetic  field  to  invert  for  (,r'vm  1)  and  (r'2,  m2). 

The  magnetic  field  is  given  by 


H  = 


R3 


3 R(R  ■  m) 

rA 


m 


(14) 


where  R  =  r  —  r'  =  x(x  -x')  +  y(y  —  y')+z(z—z')  (the  prime  coordinates  corresponding  to  the  source  and 
the  unprimed  to  the  observation)  and  m  =  xmx  +  ymy  +  zmz.  In  the  case  of  two  targets, 


h=h1+h2 


(15) 


where  H 1  and  H2  are  generated  by  the  two  targets  respectively,  and  each  obeying  Eq.  (14). 

In  the  bistatic  case,  the  dipole  m  in  Eq.  (14)  represents  the  only  source  in  the  system:  it  is  buried 
at  a  position  R  and  produces  the  magnetic  field  that  is  measured  by  the  sensor  across  the  observation 
grid.  In  the  monostatic  case,  however,  the  dipole  moment  depends  on  the  polarizability  of  the  target 
and  on  the  primary  field  at  its  location,  which  is  written  as 

m=M-Hpr,  (16) 

where  M  is  the  polarizability  tensor  and  Hpr  is  the  primary  field  at  the  target  locations.  In  the 
monostatic  configuration,  the  sensor  emits  and  receives  at  each  point  on  the  observation  grid.  At 
each  point,  the  primary  field  at  the  location  of  the  target  Hpr  needs  to  be  computed  and  is  used 
to  weight  the  polarizability  tensor  M.  This  means  that  the  dipole  moment  m  is  different  for  each 
location  of  the  sensor.  However,  the  polarizability  M  is  constant  and  needs  to  be  inverted  for,  later  to 
be  used  as  a  classification  metric  between  different  UXO. 


2.2.2  Newton  method 

In  this  section,  the  formulation  of  the  Gauss-Newton  algorithm  is  presented  for  both  monostatic  and 
bistatic  configurations,  as  well  as  when  all  the  components  of  the  magnetic  field  are  measured  or 
when  only  Hz  is  measured. 
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Bistatic  case:  Using  only  the  measured  Hz 


The  master  equation  on  which  the  iterative  scheme  is  based  is  the  standard  normal  equation 
written  as 

JJ  •  J2  •  Ax  =  JJ  •  [H?ea  -  Hczomp )  (17) 

where  x  is  the  unknown  vector  of  dimensions  [M  x  1],  Ax  =  xI+1 —  xl  is  the  unknown  and  provides 
the  updated  quantities  of  x  at  each  iteration,  Jz  is  the  Jacobian  matrix  of  dimensions  [N  x  M\,  H"lpa 
and  Hc°mp  are  the  measured  and  computed  z  components  of  the  magnetic  field,  respectively,  of 
dimensions  [Af  x  1], 

In  the  two  target  detection  case,  we  need  to  iterate  for  two  positions  and  two  momenta  so  that 


x  —  \X\,  Y\,Z\,  X2,  Y2,Z2,  Tnx±,  Ttiy\,  mzi,  mx  2,  tn  y2,niz2\ 


(18) 


where  X±  =  x  —  x'v  and  similarly  for  the  other  components.  Therefore  in  this  case,  M  =  12.  Explicitely, 
the  Jacobian  matrix  is 


h 


[Jz6l  Jz6?  Jzm^  JZfh9 


(19a) 


f  dHf 

dHf 

dH f  ) 

f  dH f 

dmx. 

dHf 

dHf 

drily. 

dHf 

dHf  \ 

dmz. 

dHf 

dX, 

dHf 

d  Yi 

dwf 

dZi 

dHf 

4,= 

d  X, 

dYi 

dZi 

>  Jzmi  ~ 

dmx. 

dmyi 

dmZi 

OH™ 

dH™ 

dH™ 

dH ™ 

dH™ 

dH ™ 

{  dX, 

dYi 

dZi  J 

{  dmH 

dmyi 

dmz.  / 

The  various  derivatives  are  given  by: 

dHz  15  XZ(R-m)  3Xmz  +  3Zmx 

~dX  ~  R7  +  R5  ’ 

dHz  _  15 YZ(Rm)  3Ymz  +  3Zmy 

~dY~  R7  +  R5  ’ 

dHz  15Z2{R-m)  6Zmz  3(R  ■  m) 

~dZ  ~  R7  +  R5  +  R5 
dHz  _  3XZ 

dmx  R5 
dHz  _  3  YZ 

dirty  R5 

dHz  _  3 Z2  1 

dm~z  ~-R5~~R3' 


(19b) 


(20a) 

(20b) 

(20c) 

(20d) 

(20e) 

(20f) 


Bistatic  case:  Using  all  components  (Hx,Hy,Hz) 


If  all  the  components  of  the  magnetic  field  are  to  be  used,  the  method  above  can  be  directly 
generalized.  Defining  Jx  and  Jy  in  the  same  way  as  in  Eq.  (19a),  the  Jacobian  and  the  right-hand 
side  term  become 


j  = 


J x 

Jy 

h 


jjj-mea  _  jzj-comp  _ 


jjj-mea  _  jfjcomp 
fjmea  fjcomp 

y  ^  y 

jjj-mea  _  fjcomp 


(21) 


where  the  vectors  over  Hx,  Hy,  Hz  denote  measurements  over  multiple  grid  points. 
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Monostatic  case  with  a  diagonal  M 


In  the  monostatic  case,  we  suppose  that  the  object  is  aligned  with  the  principal  axis  of  the  sensor 
so  that  the  polarizability  tensor  M  is  diagonal  in  the  (x,y,z)  coordinate  system.  From  Eq.  (16),  we 
therefore  write: 

m  =  xMxHpr  +  yMyHpr  +  zMzHpr.  (22) 

The  magnetic  field  is  therefore  given  by 


3 Z(R  ■  m) 
R2 


-MzHpr 


1 


3Z{XMxHpr  +  YMyHpr  +  ZMzHpr) 


R 2 


■  ~MzHpr 


(23) 


The  Jacobian  is  calculated  like  in  the  bistatic  case,  with  the  various  derivatives  given  by 

dHz  15 XZ(R-m)  3Xmz  +  3Zmx  dxmz  3Z(R  ■dxm) 

~dX  ~  Rd  +  R~5  R3~+  R5  ’ 

dHz  15  YZ(R-fh)  3Ymz  +  3Zmy  dymz  3  Z{R-dym) 

~dY=  Ri  +  R5  R^+  R~5  ’ 

dHz  15  Z2(R-?n)  6Zmz  dzmz  3  Z{R-dzm)  3  (R-m) 

~dZ~  R^  +_R5  ^3“+  ^5  +  R5 

dHz  _  3XZ 
dMx  R 5  x  ’ 

dHz  _  3YZ 

dMy  R5  y  ’ 

dHz  _  3Z2  Hpr 

dMz  R5  2  R3  ‘ 

where 

mJ1=MJ1Hpr,  T]  =  x,y,z,  (25a) 

dxm  =  xMxdxHpr  +  yMydxHpr  +  zMzdxHpr .  (25b) 

The  previous  equations  can  be  directly  used  to  build  the  Jacobian  similar  to  Eqs.  (19b),  where 
the  derivatives  with  respect  to  m  are  replaced  by  derivatives  with  respect  to  M  =  xMx  +  yMy  +  zMz. 
The  generalization  to  the  other  components  of  the  magnetic  field  is  straightforward. 


(24a) 

(24b) 

(24c) 

(24d) 

(24e) 

(240 


2.2.3  Results 

The  iterative  method  is  implemented  based  on  Eq.  (17).  The  maximum  number  of  iterations  is  set  to 
1000  (or  500,  see  below),  but  the  process  is  stopped  if  one  of  the  following  condition  applies: 

(i)  The  inverse  condition  number  of  JT  J  is  lower  than  10~10. 

(ii)  The  ten  previous  incremental  values  in  position  A R  and  in  momentum  Am  are  smaller  than 
0.1  cm  and  0.1,  respectively.  More  precisely:  at  each  iteration,  the  sum  of  increments  is  com¬ 
puted  as  dr  =  \dx\\  +  \dy\\  + . . .  +  \dz%\  and  dm  =  |dm*il  + . . .  +  \dmZ2\  and  stored  for  the  last  10 
iterations.  If  all  10  values  of  dr  and  dm  are  smaller  than  the  threshold,  it  is  decided  that  the 
algorithm  has  reached  convergence. 

The  iterative  algorithm  starts  by  making  a  random  guess  for  (i?i,mi)  and  (R^,  m2).  The  bounds 
for  the  random  intervals  of  each  components  are  as  follows:  the  lateral  positions  (x  and  y  components) 
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are  supposed  to  be  within  -50  cm  and  +50  cm,  the  depths  ( z  components)  are  supposed  to  be  within 
-20  cm  and  -50  cm,  the  momenta  are  supposed  to  be  between  -3  and  +3.  The  overall  process  is 
repeated  1000  times  (when  only  Hz  is  used)  and  500  times  (when  Hx,Hy,Hz  are  used).  The  results 
are  either  directly  averaged,  or  “bin  averaged".  In  the  latter  case,  bins  of  1  cm  (for  the  position)  or 
0.1  (for  the  momentum)  are  defined,  and  the  average  is  weighted  by  the  number  of  occurrences  in 
each  bin. 

The  algorithm  is  run  on  six  cases  corresponding  to: 

Case  1:  the  reference  case  where  two  dipoles  are  somewhat  close  to  each  other,  offset  on  the  mea¬ 
surement  grid,  and  of  different  orientations. 

Case  2:  a  variation  of  case  1  where  the  deepest  dipole  is  buried  much  deeper. 

Case  3:  a  variation  of  case  1  where  the  shallowest  dipole  has  its  momentum  drastically  reduced. 
This  could  correspond  to  a  piece  of  clutter  closer  to  the  sensor  than  the  main  object. 

Case  4:  same  as  case  1  but  the  two  objects  are  centered  on  the  measurement  grid.  The  measured 
magnetic  field  therefore  presents  less  amplitude  variations  over  the  grid. 

Case  5:  same  as  case  2  but  centered  over  the  grid. 

Case  6:  same  as  case  3  but  centered  over  the  grid. 

The  results  are  summarized  in  Tabs.  4  and  5  for  the  case  when  Hz  only  is  used,  and  in  Tabs.  6  and  7 
for  the  case  when  ( Hx,Hy,Hz )  are  used.  Results  are  also  represented  graphically  in  Figs.  7-12.  The 
following  comments  can  be  drawn: 

•  The  bin  average  usually  performs  better  than  the  direct  average. 

•  The  depth  estimation  is  rather  robust  to  noise. 

•  Paradoxically,  the  estimates  of  the  centered  cases  are  often  worse  than  those  of  the  offset  cases. 

•  It  is  not  totally  clear  that  using  the  full  magnetic  field  (Hx ,  Hy ,  Hz )  yields  better  estimates 

than  using  only  Hz.  For  example,  the  cases  where  one  object  is  deeply  buried  are  much  better 

predicted  by  using  Hz  only.  This  might  be  due  to  the  specific  dipole  moments  chosen,  and 
deserves  further  investigation. 

•  Another  poor  estimation  corresponds  to  the  situation  where  the  shallowest  object  has  a  very 
small  dipole  moment,  somewhat  analogous  to  a  cluttered  configuration.  Here  again  using  Hz 
only  yields  better  estimates  than  using  (Hx,Hy,Hz),  and  the  estimations  rapidly  degrade  with 
increasing  noise. 
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Noise 

X\ 

R 1  [cm] 

yi 

21 

X2 

R  2  [cm] 

y2 

22 

mXl 

m\ 

my  1 

mZl 

mX2 

m2 

my2 

mZ2 

Case  1 

33 

14 

-28 

13 

12 

-35 

0 

2 

1 

3 

0 

1 

0% 

32.9 

14.1 

-28.0 

13.0 

11.8 

-34.9 

0.0 

2.0 

1.0 

3.0 

0.0 

1.0 

1% 

32.9 

14.0 

-28.0 

12.9 

11.8 

-34.7 

0.0 

2.0 

1.0 

3.0 

0.0 

1.0 

2% 

32.9 

14.0 

-28.0 

12.7 

11.4 

-34.4 

0.0 

2.0 

1.0 

3.0 

-0.0 

1.0 

5% 

32.4 

14.4 

-28.6 
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11.5 

-35.6 

0.1 

2.0 

1.4 

2.5 

-0.0 

0.6 

10% 

30.5 

15.8 

-29.5 
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6.9 

-31.9 

0.6 

1.8 

2.0 

1.4 

0.0 

0.2 

Case  2 

33 

14 
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13 

12 
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0 

2 

1 

3 

0 

1 

0% 
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-28.0 

12.9 
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1.0 

3.0 
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Case  3 

33 

14 
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13 

12 
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0 

0.1 

3 

0 

1 

0% 

26.0 

13.0 

-30.3 

12.2 

13.2 

-34.6 

0.7 

-0.7 

-0.3 

2.3 

0.9 

1.4 

1% 

25.4 

13.6 

-30.0 

10.2 

11.9 

-33.7 

0.9 

-0.5 

-0.1 

2.2 

0.7 

1.2 

2% 
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14.3 

-31.0 

3.9 
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1.5 

-0.2 

0.4 
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0.4 

0.7 
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Case  4 
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2 

1 
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1 

0% 

12.6 

3.6 

-28.6 

-6.3 

1.2 

-34.1 

0.3 

1.5 

0.8 

2.7 

0.2 

1.0 

1% 

12.7 

3.6 

-28.6 

-6.3 

1.4 

-34.4 

0.3 

1.5 

0.8 

2.7 

0.2 

1.0 

2% 

12.5 

3.5 

-28.5 

-6.5 
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-34.2 
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0.2 

0.7 

Case  5 

10 

1 
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-1 

-55 

0 

2 

1 

3 

0 

1 

0% 

12.4 

3.1 

-28.4 

-6.7 

1.0 

-51.1 

0.2 

1.5 

0.8 

2.8 

0.1 

1.0 

1% 

12.1 

2.9 

-28.6 

-5.8 

0.9 

-50.1 

0.2 

1.6 

0.9 
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0.1 

0.9 

2% 

11.4 

2.3 

-28.6 

-5.8 

0.0 

-50.5 

0.2 

1.7 
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5% 
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-4.0 

-0.4 
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0.0 

0.9 

10% 

10.7 

3.4 

-29.5 

-2.8 

-0.3 
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0.4 

1.7 

1.3 

1.2 

0.1 

0.6 

Case  6 

10 

1 

-28 

-10 

-1 

-35 

0 

(hi 

3 

0 

1 

0% 

8.4 

2.9 

-28.9 

-8.1 

1.6 

-34.6 

0.5 

-0.2 

-0.0 

2.5 

0.4 

1.2 

1% 

6.1 

1.6 

-29.4 

-10.0 

0.6 

-34.5 

0.5 

-0.2 

-0.0 

2.5 

0.4 

1.1 

2% 

5.0 

2.1 

-29.5 

-8.9 

0.1 

-33.5 

0.8 

-0.0 

0.2 

2.2 

0.2 

0.9 

5% 

-1.0 

1.8 

-32.0 

-3.9 

0.2 

-24.7 

1.8 

0.1 

0.7 

1.2 

0.1 

0.4 

10% 

-2.7 

1.6 

-32.2 

-2.8 

-0.8 

-22.4 

2.0 

0.0 

0.8 

0.9 

0.1 

0.3 

Table  4:  Estimation  of  (Ri,m\)  and  (f?2, m2)  when  only  Ilz  is  used.  Values  correspond  to  direct  averages. 
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-30.7 

-10.0 

-0.7 

-34.8 

0.1 

0.1 

0.1 

2.7 

0.0 

1.0 

5% 

-8.0 

0.8 

-34.2 

-9.0 

-0.3 

-28.2 

1.8 

0.1 

0.8 

0.4 

0.0 

0.1 

10% 

-8.0 

0.8 

-34.0 

-7.4 

-0.7 

-22.3 

2.1 

0.1 

0.9 

0.2 

0.0 

0.1 

Table  5:  Estimation  of  and  (R<2,m<i)  when  only  II z  is  used.  Values  correspond  to  bin  averages. 
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Noise 

X\ 

Ri  [cm] 

yi 

z  1 

X2 

R  2  [cm] 

V2 

Z2 

mXl 

mi 

myi 

mZl 

mX2 

m2 

my2 

mZ2 

Case  1 

33 

14 

-28 

13 

12 

-35 

0 

2 

1 

3 

0 

1 

0% 

33.0 

14.0 

-28.0 

13.0 

12.0 

-35.0 

-0.0 

2.0 

1.0 

3.0 

-0.0 

1.0 

1% 

33.0 

14.1 

-28.0 

13.0 

12.1 

-34.9 

0.0 

2.0 

1.0 

3.0 

0.0 

1.0 

2% 

32.9 

14.1 

-28.0 

13.0 

11.9 

-34.9 

0.0 

2.0 

1.0 

3.0 

-0.0 

1.0 

5% 

32.6 

14.0 

-28.1 

12.7 

10.5 

-34.0 

0.1 

2.1 

1.0 

2.9 

-0.1 

1.0 

10% 

31.0 

15.1 

-29.0 

6.0 

8.1 

-33.0 

0.5 

2.0 

1.6 

2.1 

-0.0 

0.5 

Case  2 

33 

14 

-28 

13 

12 

-55 

0 

2 

1 

3 

0 

1 

0% 

33.0 

14.0 

-28.0 

13.0 

12.0 

-55.0 

0.0 

2.0 

1.0 

3.0 

-0.0 

1.0 

1% 

33.0 

14.0 

-28.0 

12.9 

11.9 

-54.9 

0.0 

2.0 

1.0 

3.0 

-0.0 

1.0 

2% 

33.0 

14.0 

-28.0 

13.1 

11.8 

-55.0 

0.0 

2.0 

1.0 

3.0 

-0.0 

1.0 

5% 

32.9 

14.0 

-28.2 

13.5 

10.2 

-52.1 

0.1 

2.0 

1.1 

2.7 

-0.1 

0.9 

10% 

32.4 

14.4 

-28.5 

7.9 

7.3 

-40.9 

0.3 

2.0 

1.3 

1.9 

-0.1 

0.7 

Case  3 

33 

14 

-28 

13 

12 

-35 

0 

(h2 

(hi 

3 

0 

1 

0% 

26.0 

12.4 

-30.7 

12.9 

13.0 

-35.4 

0.4 

-0.4 

-0.2 

2.6 

0.7 

1.4 

1% 

26.6 

13.7 

-30.1 

10.6 

10.9 

-34.3 

0.7 

-0.3 

0.0 

2.4 

0.5 

1.1 

2% 

23.8 

14.1 

-30.8 

7.2 

9.6 

-32.4 

1.1 

-0.3 

0.2 

1.9 

0.5 

0.9 

5% 

20.4 

14.2 

-32.0 

2.8 

3.5 

-25.8 

1.8 

0.0 

0.7 

1.2 

0.2 

0.5 

10% 

18.0 

13.8 

-33.3 

-0.8 

1.9 

-22.1 

2.2 

0.0 

0.8 

0.8 

0.2 

0.3 

Case  4 

10 

1 

-28 

-10 

-1 

-35 

0 

2 

1 

3 

0 

1 

0% 

12.8 

3.0 

-28.4 

-6.1 

1.3 

-35.4 

0.1 

1.6 

0.8 

2.9 

0.1 

1.1 

1% 

12.8 

3.0 

-28.4 

-6.0 

1.3 

-35.3 

0.1 

1.6 

0.8 

2.9 

0.1 

1.1 

2% 

12.7 

3.0 

-28.4 

-5.9 

1.2 

-35.3 

0.1 

1.6 

0.8 

2.9 

0.1 

1.1 

5% 

12.5 

2.9 

-28.5 

-5.8 

1.2 

-35.0 

0.2 

1.7 

0.9 

2.8 

0.0 

1.0 

10% 

11.4 

3.4 

-28.7 

-4.3 

0.1 

-32.1 

0.5 

1.7 

1.1 

2.4 

-0.0 

0.9 

Case  5 

10 

1 

-28 

-10 

-1 

-55 

0 

2 

1 

3 

0 

1 

0% 

12.1 

2.6 

-28.4 

-6.9 

0.9 

-52.7 

0.1 

1.7 

0.9 

2.9 

0.1 

1.0 

1% 

12.1 

2.6 

-28.4 

-6.7 

0.8 

-52.6 

0.1 

1.7 

0.9 

2.9 

0.1 

1.0 

2% 

11.4 

1.9 

-28.5 

-6.7 

-0.2 

-52.4 

0.1 

1.9 

1.1 

2.8 

-0.0 

0.9 

5% 

10.9 

1.8 

-28.5 

-5.9 

-0.8 

-49.4 

0.2 

1.9 

1.1 

2.6 

-0.0 

0.9 

10% 

10.6 

2.2 

-28.8 

-5.3 

0.2 

-38.8 

0.2 

1.9 

1.2 

1.9 

-0.0 

0.8 

Case  6 

10 

1 

-28 

-10 

-1 

-35 

0 

(h2 

M 

3 

0 

1 

0% 

8.0 

2.3 

-29.0 

-8.7 

0.4 

-36.5 

0.5 

0.0 

0.1 

2.4 

0.3 

1.1 

1% 

7.4 

1.0 

-30.0 

-8.1 

0.6 

-36.2 

0.3 

0.1 

0.1 

2.7 

0.3 

1.1 

2% 

6.9 

1.4 

-29.5 

-7.7 

0.9 

-35.1 

0.5 

0.1 

0.2 

2.5 

0.2 

1.0 

5% 

1.3 

2.2 

-31.3 

-5.0 

0.4 

-29.0 

1.5 

0.2 

0.8 

1.4 

0.1 

0.4 

10% 

-2.1 

1.6 

-32.9 

-1.8 

-0.7 

-22.7 

2.1 

0.3 

1.0 

0.9 

0.0 

0.2 

Table  6:  Estimation  of and  (R2, m2)  when  (Hx,Hy,Hz)  are  used.  Values  correspond  to  direct  averages. 
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Noise 

X\ 

[cm] 

yi 

21 

X2 

[cm] 

y2 

22 

mXl 

m  i 

myi 

mZl 

mx 2 

m2 

my2 

mZ2 

Case  1 

33 

14 

-28 

13 

12 

-35 

0 

2 

1 

3 

0 

1 

0% 

33.0 

14.0 

-28.0 

13.0 

12.0 

-35.0 

-0.0 

2.0 

1.0 

3.0 

-0.0 

1.0 

1% 

33.0 

14,0 

-28.0 

13.0 

12.0 

-34.9 

-0.0 

2.0 

1.0 

3.0 

0.0 

1.0 

2% 

32.9 

13.9 

-28.0 

13.1 

11.8 

-34.9 

0.0 

2.0 

1.0 

3.0 

-0.0 

1.0 

5% 

32.7 

13.9 

-28.1 

13.3 

11.8 

-34.5 

0.0 

2.0 

1.0 

2.9 

0.0 

1.0 

10% 

31.4 

14.8 

-29.3 

12.0 

11.2 

-34.1 

0.1 

1.8 

1.3 

2.0 

0.0 

0.7 

Case  2 

33 

14 

-28 

13 

12 

-55 

0 

2 

1 

3 

0 

1 

0% 

33.0 

14.0 

-28.0 

13.0 

12.0 

-55.0 

0.0 

2.0 

1.0 

3.0 

-0.0 

1.0 

1% 

33.0 

14.0 

-28.0 

12.9 

11.8 

-54.9 

0.0 

2.0 

1.0 

3.0 

-0.0 

1.0 

2% 

33.0 

14.0 

-28.0 

13.1 

12.2 

-54.8 

0.0 

2.0 

1.0 

3.0 

0.0 

1.0 

5% 

33.0 

13.9 

-28.0 

13.2 

10.9 

-53.0 

0.0 

2.0 

1.0 

2.8 

-0.0 

1.0 

10% 

32.7 

14.2 

-28.6 

11.5 

9.5 

-44.5 

0.1 

2.0 

1.1 

1.4 

-0.0 

0.5 

Case  3 

33 

14 

-28 

13 

12 

-35 

0 

0.1 

3 

0 

1 

0% 

28.3 

12.4 

-29.9 

12.9 

12.5 

-35.0 

0.2 

-0.3 

-0.2 

2.8 

0.6 

1.3 

1% 

25.7 

14.6 

-29.5 

12.8 

12.5 

-35.0 

0.2 

0.1 

0.1 

2.9 

0.1 

1.0 

2% 

19.6 

14.5 

-32.1 

12.4 

12.5 

-34.9 

0.6 

0.1 

0.3 

2.2 

0.1 

0.9 

5% 

15.8 

13.8 

-34.3 

10.6 

10.8 

-31.8 

1.9 

0.1 

0.9 

0.6 

0.0 

0.2 

10% 

14.9 

13.8 

-34.4 

6.2 

8.7 

-26.0 

2.4 

0.1 

1.0 

0.2 

0.0 

0.1 

Case  4 

10 

1 

-28 

-10 

-1 

-35 

0 

2 

1 

3 

0 

1 

0% 

10.3 

1.2 

-28.0 

-9.6 

-0.7 

-35.0 

0.0 

2.0 

1.0 

3.0 

0.0 

1.0 

1% 

10.4 

1.2 

-28.0 

-9.5 

-0.6 

-35.0 

0.0 

1.9 

1.0 

3.0 

0.0 

1.0 

2% 

10.6 

1.2 

-28.1 

-9.3 

-0.4 

-35.1 

0.0 

1.9 

0.9 

3.0 

0.0 

1.0 

5% 

11.2 

1.4 

-28.2 

-8.5 

0.1 

-34.7 

0.0 

1.7 

0.9 

2.9 

0.0 

1.1 

10% 

11.6 

2.2 

-28.9 

-6.4 

1.3 

-34.2 

0.1 

1.4 

0.8 

2.5 

0.1 

0.9 

Case  5 

10 

1 

-28 

-10 

-1 

-55 

0 

2 

1 

3 

0 

1 

0% 

10.2 

1.1 

-28.0 

-9.8 

-0.8 

-54.7 

0.0 

2.0 

1.0 

3.0 

0.0 

1.0 

1% 

10.3 

1.1 

-28.0 

-9.5 

-0.5 

-54.0 

0.0 

2.0 

1.0 

3.0 

0.0 

1.0 

2% 

10.2 

1.0 

-28.1 

-9.1 

-0.6 

-53.6 

0.0 

2.0 

1.0 

3.0 

0.0 

1.0 

5% 

10.2 

1.0 

-28.1 

-8.5 

-0.3 

-51.1 

0.0 

2.0 

1.0 

2.8 

0.0 

1.1 

10% 

9.7 

1.4 

-28.9 

-5.6 

-0.5 

-41.2 

0.1 

2.0 

1.0 

1.2 

0.0 

0.5 

Case  6 

10 

1 

-28 

-10 

-1 

-35 

0 

M 

3 

0 

1 

0% 

8.6 

1.4 

-28.3 

-10.4 

-0.4 

-34.9 

0.2 

-0.0 

0.0 

2.8 

0.2 

1.1 

1% 

8.3 

1.1 

-28.6 

-9.9 

-0.8 

-35.1 

0.0 

0.2 

0.1 

3.0 

0.0 

1.0 

2% 

6.5 

1.4 

-29.4 

-9.8 

-0.7 

-35.0 

0.0 

0.2 

0.1 

2.9 

0.0 

1.0 

5% 

-6.5 

1.1 

-33.9 

-8.9 

-0.3 

-32.6 

1.3 

0.1 

0.8 

0.7 

0.0 

0.3 

10% 

-8.2 

0.8 

-34.5 

-4.9 

0.0 

-20.9 

2.3 

0.2 

1.1 

0.1 

0.0 

0.0 

Table  7:  Estimation  of  ( Ri,m\ )  and  (R‘2,ni2)  when  (Hx,Hy,Hz)  are  used.  Values  correspond  to  bin  averages. 
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Figure  7:  Case  1:  estimations  of  and  (^2,^2)  in  the  case  of  direct  averages  and  bin  averages,  when 

both  Hz  and  ( Hx,Hy,Hz )  are  used.  The  cases  correspond  to  Tabs.  4-7.  Legend:  ‘o’  for  Ilz  straight  average;  ‘x’ 
for  Hz  bin  average;  for  ( Hx,Hy,Hz )  straight  average;  V  for  (Hx,Hy,Hz)  bin  average.  Solid  lines  correspond 
to  true  values:  x  component  is  in  blue,  y  component  in  red,  z  component  in  green. 
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(a)  Estimation  in  position. 


(b)  Estimation  in  momentum. 


Figure  8:  Case  2:  estimations  of  (Ri,mi)  and  (-R2,  m2)  in  the  case  of  direct  averages  and  bin  averages,  when 
both  Hz  and  (Hx,Hy,Hz)  are  used.  The  cases  correspond  to  Tabs.  4-7.  See  Fig.  7  for  details  on  the  legend. 
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Figure  9:  Case  3:  estimations  of  OR  i ,  m  i )  and  (R^,m^)  in  the  case  of  direct  averages  and  bin  averages,  when 
both  Hz  and  (Hx,Hy,Hz)  are  used.  The  cases  correspond  to  Tabs.  4-7.  See  Fig.  7  for  details  on  the  legend. 
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Figure  10:  Case  4:  estimations  of  CRi,mi)  and  (R?,  m2)  in  the  case  of  direct  averages  and  bin  averages,  when 
both  Hz  and  (Hx,Hy,Hz)  are  used.  The  cases  correspond  to  Tabs.  4-7.  See  Fig.  7  for  details  on  the  legend. 
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Figure  11:  Case  5:  estimations  of  (R\,m\)  and  (7? 2,  m2)  in  the  case  of  direct  averages  and  bin  averages,  when 
both  Hz  and  (Hx,Hy,Hz)  are  used.  The  cases  correspond  to  Tabs.  4-7.  See  Fig.  7  for  details  on  the  legend. 
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Figure  12:  Case  6:  estimations  of  (7?i,  m\)  and  (R 2,  m2)  in  the  case  of  direct  averages  and  bin  averages,  when 
both  Hz  and  (Hx,Hy,Hz)  are  used.  The  cases  correspond  to  Tabs.  4-7.  See  Fig.  7  for  details  on  the  legend. 
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Noise 

Case  1 

Case  2 

Case  3 

Case  4 

Case  5 

Case  6 

0% 

43.8(37.0) 

42.4  (29.0) 

38.8  (27.8) 

48.7  (38.6) 

47.6(36.0) 

43.6(38.6) 

1% 

43.2  (33.0) 

41.6(27.4) 

45.5  (38.8) 

47.3  (40.6) 

44.6  (38.8) 

39.9(36.8) 

2% 

40.3  (34.4) 

40.6  (26.0) 

47.9  (43.2) 

48.3  (41.8) 

45.4(32.6) 

43.7(35.8) 

5% 

47.0(29.6) 

54.8(30.0) 

56.0  (44.4) 

55.6(41.2) 

56.9(31.0) 

52.9(38.6) 

10% 

66.0(45.6) 

64.0  (46.0) 

58.1  (46.8) 

64.0  (50.6) 

59.1(44.2) 

55.0  (44.6) 

Table  8:  Percentage  of  divergent  cases  when  Hz  is  used  and  when  (Hx,Hy,Hz)  is  used  (the  values  correspond¬ 
ing  to  the  latter  case  are  indicated  in  parenthesis).  The  case  number  refers  to  Tabs.  4-7. 


A  statistics  in  divergence  rates  is  also  presented  in  Table  8.  It  can  be  seen  that  the  divergence 
rates  drops  sometimes  significantly  when  the  full  vectorial  magnetic  field  is  used,  as  opposed  to 
when  only  Hz  is  used.  Note  also  that  the  number  of  diverging  cases  (and  therefore  number  of  false 
detections)  depends  on  the  threshold  set  for  the  matrix  conditioning  (10“10  here).  Relaxing  this 
condition  yields  less  divergences  and  more  false  detection,  and  conversely.  However,  it  might  appear 
preferable  to  have  more  divergences  (which  can  be  discarded  in  the  decision  making  process)  for  an 
increased  probability  of  detection.  In  Table  8,  the  percentage  of  correct  detections  has  been  computed 
after  discarding  the  divergent  cases. 
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3  Conclusions  and  future  work 

Despite  the  analytical  appeal  of  the  method,  its  limitation  to  single  target  inversion  dras¬ 

tically  hinders  its  applicability  in  realistic  situations,  where  UXO  are  often  buried  alongside  clutter 
items.  In  order  to  overcome  this  limitation,  we  have  proposed  a  different  approach,  more  numer¬ 
ically  involved  but  also  of  much  broader  applicability:  the  iterative  Gauss-Newton  method.  The 
validation  of  the  method  has  been  performed  here  with  two  synthetic  targets  which  have  been  prop¬ 
erly  identified,  despite  varying  degrees  of  noise.  The  analytical  computation  of  the  Jacobian  matrices 
contributes  to  the  efficiency  of  the  method  and  make  it  a  good  candidate  for  future  investigation.  We 
therefore  suggest  additional  investigations  into  the  Gauss-Newton  algorithm  in  order  to  establish 
its  accuracy  with  real  sensor  data,  first  in  single  target  and  then  in  multi-target  configurations.  In 
particular,  the  stability  of  the  algorithm  to  the  increasing  number  of  targets  (hence  the  increasing 
number  of  unknowns)  should  be  investigated,  as  well  as  the  proper  formulation  of  the  unknown 
vector  which  could  include  a  single,  a  few,  or  all  time  channels. 
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